【题解】[HNOI2019]白兔之舞 单位根反演+MTT luoguP5293 | Qiuly's blog!

【题解】[HNOI2019]白兔之舞 单位根反演+MTT luoguP5293

Posted on | Edited on | In | Comments:

单位根反演不会啊怎么搞 $FFT$ 吧，还是了解了单位根反演后才可以搞的好吧……居然有人吐槽我说我学了 $FFT​$ 但是不会运用？！，嘤嘤嘤打击有些大……

实际上所谓的单位根反演就是这个东西：

回到题目，我们先考虑正解的简化版—— $n=1$ 的版本，我们先定义 $W=w[1][1]$ 。

现在对于每一个 $t$ 的答案显然为 $\sum_{i=0}^{L}[i\% k=t] W^i (^L_i)$

这个式子显然等于 $\sum_{i=0}^{L}[k|(i-t)] w^i (^L_i)$ 。会发现 $[k|(i-t)]$ 和上面单位根反演的 $[n|d]$ 一样，于是我们尝试将单位根反演的式子带进去。

后面的 $(\omega_k^{j} W+1)^L$ 显然可以预处理，记为 $num_j$ 。

然后发现 $-tj=\binom{t}{2}+\binom{j}{2}-\binom{t+j}{2}$

后面的式子可以用 $FFT$ 加速，但是值域太大这里需要用到 $MTT$ 。现在就有 $40$ 分了，接下来考虑 $n>1$ 的情况。

我们建矩阵，然后会发现 $n>1$ 仅会对 $num_j$ 的计算方式有变化。

我们定义一个 $begin$ 矩阵，该矩阵只有 $(0,x)$ 位置上有值且值为 $1$ ，也就是说这是白兔的起点。那么最后我们需要留下来的也就是矩阵的 $(0,y)$ ，因为只有在第二维为 $y$ 是才会计入答案。

嗯，差不多可以这样写：

1
2
3
4
5
begin.c[0][x]=1;
for(int i=0;i<k;++i) 
    num[i]=(begin*pow(w*num[i]+I,n)).c[0][y]%MOD;
/*w就是上文中的W，不过这里是矩阵*/
/*I是矩阵中的单位'1'*/

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#include <cmath>
#include <cstdio>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;

const int N=65536;
const double PI=acos(-1);   

int m,k,n,x,y,MOD,G,num[N],A[N<<2],B[N<<2];

namespace OI {
    #define F(x,i,j) for((x)=(i);(x)<=(j);++(x))
    inline ll IN() {
        char ch;bool flag=0;ll x=0;
        while(ch=getchar(),!isdigit(ch)) if(ch=='-') flag=1;
        while(isdigit(ch)) x=x*10+ch-'0',ch=getchar();
        if(flag) x=-x;return x;
    }
    struct matrix {int c[3][3];matrix(){memset(c,0,sizeof(c));}};
    matrix operator + (const matrix&a,const matrix&b) {
        matrix ans;int i,j;F(i,0,2)F(j,0,2) {
            ans.c[i][j]=a.c[i][j]+b.c[i][j];
            if(ans.c[i][j]>=MOD) ans.c[i][j]-=MOD;
        }return ans;
    }
    matrix operator * (const matrix&a,const matrix&b) {
        matrix ans;int i,j,k;F(i,0,2)F(j,0,2)F(k,0,2)
            ans.c[i][k]=(ans.c[i][k]+1ll*a.c[i][j]*b.c[j][k])%MOD;
        return ans;
    }
    matrix operator * (const matrix&a,const int&b) {
        matrix ans;int i,j;F(i,0,2)F(j,0,2)ans.c[i][j]=1ll*a.c[i][j]*b%MOD;
        return ans;
    }
    struct complex{complex(long double a=0,long double b=0){x=a,y=b;}long double x,y;};
    complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
    complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
    complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
    matrix I;
    inline int pow(int x,int y) {
        int res=1;for(;y;y>>=1,x=1ll*x*x%MOD) if(y&1) res=1ll*res*x%MOD;
        return res%MOD;
    }
    inline matrix pow(matrix x,int y) {
        matrix res=I;for(;y;y>>=1,x=x*x) if(y&1) res=res*x;
        return res;
    }
}using namespace OI;

namespace MTT {
    #define BLOCK 32768
    int limit=1,cnt=0,filp[N<<2],Ans[N<<2];
    complex A1[N<<2],B1[N<<2],A2[N<<2],B2[N<<2],X[N<<2],omg[N<<2];
    inline void fft(complex *f,short inv){
        for(int i=0;i<limit;++i)if(i<filp[i])std::swap(f[i],f[filp[i]]);
        for(int p=1;p<limit;p<<=1)
            for(complex *a=f;a!=f+limit;a+=(p<<1))
                for(int l=0;l<p;++l){
                    complex t=a[l+p]*omg[limit/(p<<1)*l];
                    a[l+p]=a[l]-t,a[l]=a[l]+t;
                }			
    }
    inline void mtt(int *A,int *B){
        while(limit<(k*3+5)) limit<<=1,++cnt;
        for(int i=0;i<limit;++i) filp[i]=(filp[i>>1]>>1)|((i&1)<<(cnt-1));
        for(int i=0;i<limit;++i) A1[i].x=A[i]&(BLOCK-1),A2[i].x=A[i]>>15;
        for(int i=0;i<limit;++i) B1[i].x=B[i]&(BLOCK-1),B2[i].x=B[i]>>15;
        for(int i=0;i<limit;++i) omg[i]=(complex){cos(i*PI*2/limit),sin(i*PI*2/limit)};
        fft(A1,1),fft(B1,1);fft(A2,1),fft(B2,1);
        for(int i=0;i<limit;++i){
            complex a1=A1[i],a2=A2[i],b1=B1[i],b2=B2[i];
            A1[i]=a1*b1,A2[i]=a2*b2,B1[i]=a1*b2,B2[i]=a2*b1;
        }
        for(int i=0;i<limit;++i) omg[i]=(complex){cos(i*PI*2/limit),-sin(i*PI*2/limit)};
        fft(A1,-1),fft(B1,-1);fft(A2,-1),fft(B2,-1);
        for(int i=0;i<limit;++i)
            A1[i].x/=limit,A2[i].x/=limit,B1[i].x/=limit,B2[i].x/=limit;
        for(int i=0;i<limit;++i)
            Ans[i]=((ll)(A1[i].x+0.5)%MOD+1073741824ll*((ll)(A2[i].x+0.5)%MOD)%MOD+
            32768ll*((ll)(B1[i].x+0.5)%MOD)%MOD+32768ll*((ll)(B2[i].x+0.5)%MOD)%MOD)%MOD;
    }
}using namespace MTT;

int divisor[105],tot;
inline int get_G() {/*获取原根*/
    for(int i=2,S=MOD-1;i<=S;++i) 
        if(S%i==0) {divisor[++tot]=i;while(!(S%i)) S/=i;}
    for(int g=2;;++g) {
        bool ok=true;
        for(int j=1;j<=tot;++j) if(pow(g,(MOD-1)/divisor[j])==1) {ok=false;break;}
        if(ok) return g;
    }
}

matrix w,s;

int main() {
    I.c[0][0]=I.c[1][1]=I.c[2][2]=1;
    m=IN(),k=IN(),n=IN(),x=IN(),y=IN(),MOD=IN();--x,--y;
    /*num其实就是上文中的单位根，这里预处理一下计算方便些*/
    num[0]=1,num[1]=pow(G=get_G(),(MOD-1)/k);
    for(int i=2;i<k;++i) num[i]=1ll*num[1]*num[i-1]%MOD;
    for(int i=0;i<m;++i) for(int j=0;j<m;++j) w.c[i][j]=IN();
    for(int i=0;i<(k<<1|1);++i) A[i]=num[(k-1ll*i*(i-1)/2%k)%k];
    s.c[0][x]=1;
    for(int i=0;i<k;++i) B[i]=1ll*num[1ll*i*(i-1)/2%k]*(s*pow(w*num[i]+I,n)).c[0][y]%MOD;
    /*计算后面两个多项式的值*/
    reverse(B,B+k+1),mtt(A,B);
    int invk=pow(k,MOD-2);
    for(int i=0;i<k;++i) printf("%lld\n",1ll*Ans[i+k]*invk%MOD*num[1ll*i*(i-1)/2%k]%MOD);
    /*计算答案*/
    return 0;
}

本文标题:【题解】[HNOI2019]白兔之舞 单位根反演+MTT luoguP5293

文章作者:Qiuly

发布时间:2019年04月17日 - 00:00

最后更新:2019年04月18日 - 09:53

原始链接:http://qiulyblog.github.io/2019/04/17/[题解]luoguP5293/

许可协议: 署名-非商业性使用-禁止演绎 4.0 国际 转载请保留原文链接及作者。